# ifndef _BSPLINE_H_
# define _BSPLINE_H_

# include"InterpolationStruct.h"
# include"Bbase.h"

/**
* @brief 定义根据节点和类型，边值条件实现 B 样条插值拟合的求解器类
* 
*/
class Bspline
{
private:
    int degree = 3;                   // 次数
    int type = S32;                   // 使用的插值类型，默认为 S32
    std::vector<struct point> points; // 存储节点信息
    Bbase basis;                      // 存储 B 样条基函数
    std::vector<double> coeff;        // 存储系数

    // S32 的边值条件
    int boundary = Nature;
    double left = 0;
    double right = 0;

    void BsplineS01();                 // S01 构造函数
    void BsplineS32Complete();         // 三次函数分段插值 S32 Complete 类型
    void BsplineS32Specified_Nature(); // 三次函数分段插值 S32 Specified 和 Nature 类型
    void BsplineS32Periodic();         // 三次函数分段插值 S32 Periodic 类型

    void BsplineCardinalS21();         // Cardinal S21 类型插值
    void BsplineCardinalS32();         // Cardinal S32 类型插值

    double getMul(const double & _a, const double & _b, const double & _c) const; // 获取乘数
    void TypeForm();        // 根据给定参数进行插值结构构造
    bool judgeXray() const; // 判断输入数据格式是否正确
    int search_x(const double & _x) const; // 寻找 x 所在的区间段

public:
    // 初始化构造
    Bspline();
    // 按照节点信息构造
    Bspline(const int & _type, const std::vector<struct point> & _points, 
           const int & _boundary = Nature, const double & _left = 0, const double & _right = 0);
    // 按照节点向量信息构造
    Bspline(const int & _type, const std::vector<double> & x, const std::vector<double> & y,
           const int & _boundary = Nature, const double & _left = 0, const double & _right = 0);
    // 按照自变量值以及给定拟合函数进行构造
    Bspline(const int & _type, const std::vector<double> & x, const Function & func,
           const int & _boundary = Nature);

    // 重载操作符 =
    void operator= (Bspline && B);
    void operator= (Bspline & B);

    double solve(const double & _x) const;       // 根据给定 x 值计算插值多项式计算值
    double diff(const double & _x) const;        // 根据给定 x 值计算插值多项式导数值
    double second_diff(const double & _x) const; // 根据给定 x 值计算插值多项式二阶导数值

    void PrintCoff() const;               // 打印系数
    void PrintPoint() const;              // 打印节点信息
};

// 判断是否有输入错的数据
bool Bspline::judgeXray() const
{
    bool judge = true;
    if (points.size() == 0)
    {
        judge = false;
        return judge;
    }

    for (int i = 0; i < points.size() - 1; i++)
    {
        if (points[i].x >= points[i+1].x)
        {
            judge = false;
            break;
        }
    }
    return judge;
};

// 根据给定类型构造对应的求解器
void Bspline::TypeForm()
{
    // 先对数组进行排序
    std::sort(points.begin(), points.end(), xUpsort);

    // 判断输入点是否横坐标有重复
    if(!judgeXray())
    {
        std::cout << "Values of x are wrong." << std::endl;
        return;
    }

    // 初始化基函数
    std::vector<double> xpoint;
    for(int i = 0; i < points.size(); i++)
        xpoint.push_back(points[i].x);
    if(type != S01 && type != S32 && type != S21Car && type != S32Car)
    {
        std::cout << "The type can not be formed." << std::endl;
        return;
    }
    
    switch (type)
    {
    case S01: degree = 1; break;
    case S21Car: degree = 2; break;
    case S32: case S32Car: degree = 3; break;
    default: break;
    }
    Bbase B(degree, xpoint);
    basis = B;

    // 根据类型选择构造系数矩阵
    switch (type)
    {
    case S01: BsplineS01(); break;
    case S32: 
        switch(boundary)
        {
            case Complete: BsplineS32Complete(); break;
            case Specified: case Nature: BsplineS32Specified_Nature(); break;
            case Periodic: BsplineS32Periodic(); break;
            default:
                std::cout << "The boundary condition is false." << std::endl;
                break;
        } break;
    case S21Car: BsplineCardinalS21(); break;
    case S32Car: BsplineCardinalS32(); break;
    default: 
        break;
    }
};

// 初始化构造
Bspline::Bspline()
{
    struct point p1{0,0}, p2{1,0};
    points.push_back(p1);
    points.push_back(p2);
    TypeForm();
}

// 按照节点信息构造
Bspline::Bspline(const int & _type, const std::vector<struct point> & _points, 
               const int & _boundary, const double & _left, const double & _right) :
    type{_type}, points{_points}, boundary{_boundary}, left{_left}, right{_right}
{
    TypeForm();
};

// 按照向量信息构造
Bspline::Bspline(const int & _type, const std::vector<double> & x, const std::vector<double> & y,
               const int & _boundary, const double & _left, const double & _right) :
    type{_type}, boundary{_boundary}, left{_left}, right{_right}
{
    // 判断点列是否数据对齐
    if (x.size() != y.size())
    {
        std::cout << "Information of points is not completed." << std::endl;
        return;
    }
    for (int i = 0; i < x.size(); i++)
    {
        struct point node{x[i], y[i]};
        points.push_back(node);
    }
    
    TypeForm();
};

// 根据节点与函数进行构造
Bspline::Bspline(const int & _type, const std::vector<double> & x, const Function & func,
               const int & _boundary) :
    type{_type}, boundary{_boundary}
{
    // Compute the boundary condition of Special case
    if(type == S21Car)
    {
        left = func(x[0] - 0.5);
        right = func(x[x.size() - 1] + 0.5);
        struct point node0{x[0]-0.5, left};
        points.push_back(node0);
        for (int i = 0; i < x.size(); i++)
        {
            struct point node{x[i] + 0.5, func(x[i])};
            points.push_back(node);
        }
    }
    else if(boundary == Complete || type == S32Car)
    {
        left = func.diff(x[0]);
        right = func.diff(x[x.size() - 1]);
    }
    else if(boundary == Specified)
    {
        left = func.second_diff(x[0]);
        right = func.second_diff(x[x.size() - 1]);
    }

    if(type != S21Car)
    {
        for (int i = 0; i < x.size(); i++)
        {
            struct point node{x[i], func(x[i])};
            points.push_back(node);
        }
    }
    
    TypeForm();
};

// 重载操作符 =
void Bspline::operator=(Bspline && B)
{
    degree = B.degree;
    type = B.type;

    std::vector<struct point> _points(B.points);
    points.swap(_points);
    basis = B.basis;
    std::vector<double> _coeff(B.coeff);
    coeff.swap(_coeff);

    // S32 的边值条件
    boundary = B.boundary;
    left = B.left;
    right = B.right;
}

void Bspline::operator=(Bspline & B)
{
    degree = B.degree;
    type = B.type;

    points.assign(B.points.begin(), B.points.end());
    basis = B.basis;
    coeff.assign(B.coeff.begin(), B.coeff.end());

    // S32 的边值条件
    boundary = B.boundary;
    left = B.left;
    right = B.right;
}


// 获取乘数
double Bspline::getMul(const double & _a, const double & _b, const double & _c) const
{
    double a = fabs(_a), b = fabs(_b), c = fabs(_c);
    double mul = 1.0;
    if(a > 0 && a < 1.0) mul = 1.0 / a;
    if(b > 0 && b < 1.0 / mul) mul = 1.0 / b;
    if(c > 0 && c < 1.0 / mul) mul = 1.0 / c;
    return mul;
}

// 线性分段插值
void Bspline::BsplineS01()
{
    for(int i = 0; i < points.size(); i++)
        coeff.push_back(points[i].y);
};

// 三次函数分段插值 S32 Complete 类型
void Bspline::BsplineS32Complete()
{
    int n = points.size();
    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n+2, n+2);
    Eigen::VectorXd b(n+2);
    for (int i = 0; i < n + 2; i++)
    {
        double mul = 1.0;
        if(i == 0)
        {
            A(0,0) = basis.diff(-2, points[0].x); A(0,1) = basis.diff(-1, points[0].x);
            A(0,2) = basis.diff(0, points[0].x); b(0) = left;
            mul = getMul(A(0,0), A(0,1), A(0,2));
            A(0,0) *= mul; A(0,1) *= mul; A(0,2) *= mul; b(0) *= mul;
        }
        else if(i == n + 1)
        {
            A(n+1, n-1) = basis.diff(n-3, points[n-1].x); A(n+1, n) = basis.diff(n-2, points[n-1].x);
            A(n+1, n+1) = basis.diff(n-1, points[n-1].x); b(n+1) = right;
            mul = getMul(A(n+1, n-1), A(n+1, n), A(n+1, n+1));
            A(n+1, n-1) *= mul; A(n+1, n) *= mul; A(n+1, n+1) *= mul; b(n+1) *= mul;
        }
        else
        {
            A(i, i-1) = basis(i-3, points[i-1].x); A(i,i) = basis(i-2, points[i-1].x);
            A(i, i+1) = basis(i-1, points[i-1].x); b(i) = points[i-1].y;
            mul = getMul(A(i, i-1), A(i, i), A(i, i+1));
            A(i, i-1) *= mul; A(i, i) *= mul; A(i, i+1) *= mul; b(i) *= mul;
        }
    }
    
    Eigen::VectorXd m = A.colPivHouseholderQr().solve(b);
    for (int i = 0; i < n+2; i++)
        coeff.push_back(m(i));
};

// 三次函数分段插值 S32 Specified 和 Nature 类型
void Bspline::BsplineS32Specified_Nature()
{
    int n = points.size();
    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n+2, n+2);
    Eigen::VectorXd b(n+2);
    for (int i = 0; i < n + 2; i++)
    {
        double mul = 1.0;
        if(i == 0)
        {
            A(0,0) = basis.second_diff(-2, points[0].x); A(0,1) = basis.second_diff(-1, points[0].x);
            A(0,2) = basis.second_diff(0, points[0].x); b(0) = left;
            mul = getMul(A(0,0), A(0,1), A(0,2));
            A(0,0) *= mul; A(0,1) *= mul; A(0,2) *= mul; b(0) *= mul;
        }
        else if(i == n + 1)
        {
            A(n+1, n-1) = basis.second_diff(n-3, points[n-1].x); A(n+1, n) = basis.second_diff(n-2, points[n-1].x);
            A(n+1, n+1) = basis.second_diff(n-1, points[n-1].x); b(n+1) = right;
            mul = getMul(A(n+1, n-1), A(n+1, n), A(n+1, n+1));
            A(n+1, n-1) *= mul; A(n+1, n) *= mul; A(n+1, n+1) *= mul; b(n+1) *= mul;
        }
        else
        {
            A(i, i-1) = basis(i-3, points[i-1].x); A(i,i) = basis(i-2, points[i-1].x);
            A(i, i+1) = basis(i-1, points[i-1].x); b(i) = points[i-1].y;
            mul = getMul(A(i, i-1), A(i, i), A(i, i+1));
            A(i, i-1) *= mul; A(i, i) *= mul; A(i, i+1) *= mul; b(i) *= mul;
        }
    }
    
    Eigen::VectorXd m = A.colPivHouseholderQr().solve(b);

    for (int i = 0; i < n+2; i++)
        coeff.push_back(m(i));
};

// 三次函数分段插值 S32 Periodic 类型
void Bspline::BsplineS32Periodic()
{
    // Judge the periodic conditon.
    if(fabs(points[0].y - points[points.size() - 1].y) > 0.000001)
    {
        std::cout << "Cannot satisfy the periodic condition." << std::endl;
        return;
    }
    points[points.size() - 1].y = points[0].y;

    // Calculate
    int n = points.size();
    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n+2, n+2);
    Eigen::VectorXd b(n+2);
    for (int i = 0; i < n + 2; i++)
    {
        double mul = 1.0;
        if(i == 0)
        {
            A(0,0) = basis.diff(-2, points[0].x); A(0,1) = basis.diff(-1, points[0].x);
            A(0,2) = basis.diff(0, points[0].x);  A(0,n-1) = - basis.diff(n-3, points[n-1].x);
            A(0,n) = - basis.diff(n-2, points[n-1].x); A(0,n+1) = - basis.diff(n-1, points[n-1].x);
            b(0) = 0;
            mul = getMul(A(0,0), A(0,1), A(0,2));
            double mul1 = getMul(A(0,n-1), A(0,n), A(0,n+1));
            if(mul1 > mul) mul = mul1;
            A(0,0) *= mul; A(0,1) *= mul; A(0,2) *= mul;
            A(0,n-1) *= mul; A(0,n) *= mul; A(0,n+1) *= mul;
        }
        else if(i == n + 1)
        {
            A(n+1,0) = basis.second_diff(-2, points[0].x); A(n+1,1) = basis.second_diff(-1, points[0].x);
            A(n+1,2) = basis.second_diff(0, points[0].x);  A(n+1,n-1) = - basis.second_diff(n-3, points[n-1].x);
            A(n+1,n) = - basis.second_diff(n-2, points[n-1].x); A(n+1,n+1) = - basis.second_diff(n-1, points[n-1].x);
            b(n+1) = 0;
            mul = getMul(A(n+1,0), A(n+1,1), A(n+1,2));
            double mul1 = getMul(A(n+1,n-1), A(n+1,n), A(n+1,n+1));
            if(mul1 > mul) mul = mul1;
            A(n+1,0) *= mul; A(n+1,1) *= mul; A(n+1,2) *= mul;
            A(n+1,n-1) *= mul; A(n+1,n) *= mul; A(n+1,n+1) *= mul;
        }
        else
        {
            A(i, i-1) = basis(i-3, points[i-1].x); A(i,i) = basis(i-2, points[i-1].x);
            A(i, i+1) = basis(i-1, points[i-1].x); b(i) = points[i-1].y;
            mul = getMul(A(i, i-1), A(i, i), A(i, i+1));
            A(i, i-1) *= mul; A(i, i) *= mul; A(i, i+1) *= mul; b(i) *= mul;
        }
    }
    
    Eigen::VectorXd m = A.colPivHouseholderQr().solve(b);

    for (int i = 0; i < n+2; i++)
        coeff.push_back(m(i));
};

// Cardinal S21 类型插值
void Bspline::BsplineCardinalS21()
{
    int n = points.size() - 1;
    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n, n);
    Eigen::VectorXd b(n);
    for (int i = 0; i < n; i++)
    {
        if(i == 0)
        {
            A(0,0) = 5; A(0,1) = 1; b(0) = 8 * points[1].y - 2 * left;
        }
        else if(i == n - 1)
        {
            A(n-1, n-2) = 1; A(n-1, n-1) = 5;
            b(n-1) = 8 * points[n].y - 2 * right;
        }
        else
        {
            A(i, i-1) = 1; A(i, i) = 6; A(i, i+1) = 1;
            b(i) = 8 * points[i + 1].y;
        }
    }
    Eigen::VectorXd M = A.colPivHouseholderQr().solve(b);

    double a0 = 2 * left - M(0), an = 2 * right - M(n-1);
    coeff.push_back(a0);
    for (int i = 0; i < n; i++)
        coeff.push_back(M(i));
    coeff.push_back(an);
};

// Cardinal S32 类型插值
void Bspline::BsplineCardinalS32()
{
    int n = points.size();
    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n, n);
    Eigen::VectorXd b(n);
    for (int i = 0; i < n; i++)
    {
        if(i == 0)
        {
            A(0,0) = 2; A(0,1) = 1; b(0) = 3 * points[0].y + left;
        }
        else if(i == n - 1)
        {
            A(n-1, n-2) = 1; A(n-1, n-1) = 2;
            b(n-1) = 3 * points[n-1].y - right;
        }
        else
        {
            A(i, i-1) = 1; A(i, i) = 4; A(i, i+1) = 1;
            b(i) = 6 * points[i].y;
        }
    }
    Eigen::VectorXd M = A.colPivHouseholderQr().solve(b);

    double a0 = M(1) - 2 * left, an = M(n-2) + 2 * right;
    coeff.push_back(a0);
    for (int i = 0; i < n; i++)
        coeff.push_back(M(i));
    coeff.push_back(an);
};

// 打印节点信息
void Bspline::PrintPoint() const
{
    for (int i = 0; i < points.size(); i++)
    {
        std::cout << points[i].x << "," << points[i].y << " ";
    }
    std::cout << std::endl;
};

// 打印系数信息
void Bspline::PrintCoff() const
{
    PrintVec(coeff);
};


// 寻找 x 所在的插值区间 左闭右开，右端点单独取出来计算
int Bspline::search_x(const double & _x) const
{
    int le = 0, ri = points.size() - 1, index = 0;
    if(_x >= points[ri].x)
        index = ri;
    else
    {
        while (le < ri)
        {
            int mid = (le + ri) / 2;
            if (_x == points[mid].x)
            {
                index = mid; break;
            }
            else if (_x < points[mid].x) ri = mid;
            else
            {
                if(_x <= points[mid + 1].x)
                {
                    index = mid; break;
                }
                else le = mid + 1;
            }
        }
        if (le == ri) index = le;
    }

    return index;
};


// 根据给定 x 值计算插值函数值
double Bspline::solve(const double & _x) const
{
    if(_x < points[0].x || _x > points[points.size() - 1].x + 0.00001)
    {
        std::cout << "Out of range." << std::endl;
        return 0;
    }

    int index = search_x(_x);
    double output = 0;
    if(degree == 1)
        output = coeff[index] * basis(index, _x) + coeff[index + 1] * basis(index + 1, _x);
    else if(degree == 2)
    {
        output = coeff[index] * basis(index - 1, _x) + coeff[index + 1] * basis(index, _x)
                 + coeff[index + 2] * basis(index + 1, _x);
    }
    else
        output = coeff[index] * basis(index - 2, _x) + coeff[index + 1] * basis(index - 1, _x)
                 + coeff[index + 2] * basis(index, _x) + coeff[index + 3] * basis(index + 1, _x);

    return output;
};

// 根据给定 x 值计算插值函数导数值
double Bspline::diff(const double & _x) const
{
    if(_x < points[0].x || _x > points[points.size() - 1].x + 0.00001)
    {
        std::cout << "Out of range." << std::endl;
        return 0;
    }

    int index = search_x(_x);
    double output = 0;
    if(degree == 1)
        output = coeff[index] * basis.diff(index, _x) + coeff[index + 1] * basis.diff(index + 1, _x);
    else if(degree == 2)
    {
        output = coeff[index] * basis.diff(index - 1, _x) + coeff[index + 1] * basis.diff(index, _x)
                 + coeff[index + 2] * basis.diff(index + 1, _x);
    }
    else
        output = coeff[index] * basis.diff(index - 2, _x) + coeff[index + 1] * basis.diff(index - 1, _x)
                 + coeff[index + 2] * basis.diff(index, _x) + coeff[index + 3] * basis.diff(index + 1, _x);

    return output;
};

// 根据给定 x 值计算插值函数二阶导数值
double Bspline::second_diff(const double & _x) const
{
    if(_x < points[0].x || _x > points[points.size() - 1].x + 0.00001)
    {
        std::cout << "Out of range." << std::endl;
        return 0;
    }

    double output = 0;
    int index = search_x(_x);
    if(degree == 2)
    {
        output = coeff[index] * basis.second_diff(index - 1, _x) + coeff[index + 1] * basis.second_diff(index, _x)
                 + coeff[index + 2] * basis.second_diff(index + 1, _x);
    }
    else if(degree == 3)
    {
        output = coeff[index] * basis.second_diff(index - 2, _x) + coeff[index + 1] * basis.second_diff(index - 1, _x)
                 + coeff[index + 2] * basis.second_diff(index, _x) + coeff[index + 3] * basis.second_diff(index + 1, _x);
    }
        
    return output;
};

# else
// Do nothing.
# endif